---
title: "Supporting Information Analysis for Assessing Public Value Failure in Government Adoption of AI"
author: "Daniel Schiff, Kaylyn Jackson Schiff, and Patrick Pierson"
date: "2021"
output: pdf_document
editor_options: 
  chunk_output_type: console
---

#####Setup Chunk#####
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(estimatr)
library(dplyr)
library(qwraps2)
library(kableExtra)
library(stringr)
library(xtable)
library(ggpubr)
library(tidyr)
library(forcats)
library(stargazer)
library(sandwich)
library(dotwhisker)
library(rmarkdown)
```

#####Read in the Clean Data####
```{r read in data}
load("Data/ads_clean.RData")
```

#####Tests for Wing Order Effects####
```{r test for wing order effects}
#Assess significance of wing order at 5% level
wingorder1 <- lm(Feeling ~ Wing_Order, data = ads)
wingorder2 <- lm(Trust ~ Wing_Order, data = ads)
wingorder3 <- lm(Quality ~ Wing_Order, data = ads)
wingorder4 <- lm(Impact ~ Wing_Order, data = ads)

f_stat_p_value = function(x){
  1 - as.numeric(pf(summary(x)$fstatistic[1], summary(x)$fstatistic[2], summary(x)$fstatistic[3]))
}

f.wingorder1 <- f_stat_p_value(wingorder1)
f.wingorder2 <- f_stat_p_value(wingorder2)
f.wingorder3 <- f_stat_p_value(wingorder3)
f.wingorder4 <- f_stat_p_value(wingorder4)

f.wingorder.frame <- rbind(f.wingorder1, f.wingorder2, f.wingorder3, f.wingorder4)
row.names(f.wingorder.frame) <- c("Feeling", "Trust", "Quality", "Impact")
colnames(f.wingorder.frame) <- c("p-value")

#Results of f-tests
f.wingorder.frame
```

####Reproduce Figure SI1####
```{r figure SI1}
#Prepare simulation
alpha <- 0.05 # Standard significance level 
possible.n <- seq(from=500, to=2000, by=50) # The sample sizes we'll be considering 
sims <- 1000 # Number of simulated treatment allocations to conduct for each N 

powers <- rep(NA, length(possible.n)) # Empty object to collect lm simulation estimates 

#Loop to vary the number of subjects
for (j in 1:length(possible.n)){ N <- possible.n[j] # Pick the sample size
set.seed(30308)
Y0 <- rnorm(n=N, mean=0, sd=2) # control potential outcome 
tau1 <- 0.5 # Hypothesize treatment effect for Fairness treatment group
Y1 <- Y0 - tau1 # treatment potential outcome for Fairness treatment group
tau2 <- 0.5 # Hypothesize treatment effect for Transparency treatment group
Y2 <- Y0 - tau2 # treatment potential outcome for Transparency treatment group
tau3 <- 0.5 # Hypothesize treatment effect for Responsiveness treatment group
Y3 <- Y0 - tau3 # treatment potential outcome for Responsiveness treatment group

significant.lm <- rep(NA, sims) # Empty object to count significant experiments 

#Inner loop to conduct experiments "sims" times over for each N
for (i in 1:sims){
  Z.sim <-  sample(0:3, N, replace = T, prob = rep(1/4,4)) # Do a random assignment 
  Y.sim <- case_when(Z.sim == 0 ~ Y0, Z.sim == 1 ~ Y1, Z.sim == 2 ~ Y2, Z.sim == 3 ~ Y3) 
  # Reveal outcomes according to assignment 
  
  fit.lm <- lm_robust(Y.sim ~ factor(Z.sim)) # Do analysis (Simple regression)
  
  p.lm1 <- summary(fit.lm)$coefficients[2,4] # Extract p-values  
  p.lm2 <- summary(fit.lm)$coefficients[3,4] # Extract p-values 
  p.lm3 <- summary(fit.lm)$coefficients[4,4] # Extract p-values 
  significant.lm[i] <- (p.lm1 <= alpha & p.lm2 <= alpha & p.lm3 <= alpha) # Determine significance according to p <= 0.05
}
powers[j] <- mean(significant.lm) # store average success rate (power) for each N 
}

#Create figure
png("figures/fig_SI1.png")
plot(possible.n, powers, type = "b", main = "Power by Sample Size: Primary Hypothesis", 
     xlab = "Sample Size", ylab = "Power")
points(possible.n, powers, 
       pch = 19, col = "gray", bg = "gray", cex = 1.3)
abline(h = .8, lty = 2)
dev.off()
```

####Reproduce Figure SI2####
```{r figure SI2}
#Prepare simulation
#Loop to vary the number of subjects
for (j in 1:length(possible.n)){ N <- possible.n[j] # Pick the sample size
set.seed(30308)
Y0 <- rnorm(n=N, mean=0, sd=2) # control potential outcome 
tau1 <- 0.5 # Hypothesize treatment effect
Y1 <- Y0 - tau1 # treatment potential outcome
tau2 <- 0.5 # Hypothesize treatment effect
Y2 <- Y0 - tau2 # treatment potential outcome
tau3 <- 0.5 # Hypothesize treatment effect
Y3 <- Y0 - tau3 # treatment potential outcome

tau4 <- 0 # Hypothesize treatment effect
Y4 <- Y0 - tau4 # treatment potential outcome
tau5 <- 0.5 # Hypothesize treatment effect
Y5 <- Y0 - tau5 # treatment potential outcome
tau6 <- 0.5 # Hypothesize treatment effect
Y6 <- Y0 - tau6 # treatment potential outcome
tau7 <- 0.5 # Hypothesize treatment effect
Y7 <- Y0 - tau7 # treatment potential outcome

significant.lm <- rep(NA, sims) # Empty object to count significant experiments 

#Inner loop to conduct experiments "sims" times over for each N
for (i in 1:sims){
  Z.sim <-  sample(0:7, N, replace = T, prob = rep(1/8,8)) # Do a random assignment 
  Y.sim <- case_when(Z.sim == 0 ~ Y0, Z.sim == 1 ~ Y1, Z.sim == 2 ~ Y2, Z.sim == 3 ~ Y3,
                     Z.sim == 4 ~ Y4, Z.sim == 5 ~ Y5, Z.sim == 6 ~ Y6, Z.sim == 7 ~ Y7) 
  # Reveal outcomes according to assignment 

  fit.lm <- lm_robust(Y.sim ~ factor(Z.sim)) # Do analysis (Simple regression)
  
  p.lm1 <- summary(fit.lm)$coefficients[2,4] # Extract p-values  
  p.lm2 <- summary(fit.lm)$coefficients[3,4] # Extract p-values  
  p.lm3 <- summary(fit.lm)$coefficients[4,4] # Extract p-values  
  p.lm4 <- summary(fit.lm)$coefficients[5,4] # Extract p-values
  p.lm5 <- summary(fit.lm)$coefficients[6,4] # Extract p-values  
  p.lm6 <- summary(fit.lm)$coefficients[7,4] # Extract p-values  
  p.lm7 <- summary(fit.lm)$coefficients[8,4] # Extract p-values  
 
  significant.lm[i] <- (p.lm1 <= alpha | p.lm2 <= alpha | p.lm3 <= alpha | p.lm4 <= alpha |
                          p.lm5 <= alpha | p.lm6 <= alpha | p.lm7 <= alpha) 
  # Determine significance according to p <= 0.05
}
  powers[j] <- mean(significant.lm) # store average success rate (power) for each N 
}

#Create figure
png("figures/fig_SI2.png")
plot(possible.n, powers, type = "b", main = "Power by Sample Size: Policy Sector Hypothesis", 
     xlab = "Sample Size", ylab = "Power")
points(possible.n, powers, 
       pch = 19, col = "gray", bg = "gray", cex = 1.3)
abline(h = .8, lty = 2)
dev.off()
```

#####Reproduce Table SI1####
```{r table SI1}
balance_table <-
  list("Race" =
         list("White" = ~ qwraps2::mean_sd(.data$Race == "White", denote_sd="paren"),
              "Black"  = ~ qwraps2::mean_sd(.data$Race == "Black or African American", denote_sd="paren"),
              "Other"  = ~ qwraps2::mean_sd(.data$Race == "Other Race", denote_sd="paren")),
       "Gender" =
         list("Male" = ~ qwraps2::mean_sd(.data$Gender == "Male", denote_sd="paren"),
              "Female"  = ~ qwraps2::mean_sd(.data$Gender == "Female", denote_sd="paren"),
              "Other"  = ~ qwraps2::mean_sd(.data$Gender == "Other", denote_sd="paren")),
       "Educ" = 
         list("High school or less" = ~ qwraps2::mean_sd(.data$Educ == "High school graduate or less", denote_sd="paren"),
              "Some college" = ~ qwraps2::mean_sd(.data$Educ == "Some college", denote_sd="paren"),
              "Bachelor's or higher" = ~ qwraps2::mean_sd(.data$Educ == "Bachelor's or higher", denote_sd="paren")),
       "Age" = 
         list("Gen Z" = ~ qwraps2::mean_sd(.data$Age == "Gen Z", denote_sd="paren"),
              "Millenial" = ~ qwraps2::mean_sd(.data$Age == "Millenial", denote_sd="paren"),
              "Gen X" = ~ qwraps2::mean_sd(.data$Age == "Gen X", denote_sd="paren"),
              "Boomers" = ~ qwraps2::mean_sd(.data$Age == "Boomers", denote_sd="paren"),
              "Silent" = ~ qwraps2::mean_sd(.data$Age == "Silent", denote_sd="paren")),
       "Pol" = 
         list("Independent" = ~ qwraps2::mean_sd(.data$Pol == "Independent", denote_sd="paren"),
              "Liberal" = ~ qwraps2::mean_sd(.data$Pol == "Liberal", denote_sd="paren"),
              "Conservative" = ~ qwraps2::mean_sd(.data$Pol == "Conservative", denote_sd="paren")),
       "Knowledge" = 
         list("A lot" = ~ qwraps2::mean_sd(.data$Knowledge == "A lot", denote_sd="paren"),
              "A moderate amount" = ~ qwraps2::mean_sd(.data$Knowledge == "A moderate amount", denote_sd="paren"),
              "A little" = ~ qwraps2::mean_sd(.data$Knowledge == "A little", denote_sd="paren"),
              "None at all" = ~ qwraps2::mean_sd(.data$Knowledge == "Not at all", denote_sd="paren")),
       "Marital Status" = 
         list("Married" = ~ qwraps2::mean_sd(.data$Marital_Status == "Married", denote_sd="paren"),
              "Widowed" = ~ qwraps2::mean_sd(.data$Marital_Status == "Widowed", denote_sd="paren"),
              "Divorced" = ~ qwraps2::mean_sd(.data$Marital_Status == "Divorced", denote_sd="paren"),
              "Separated" = ~ qwraps2::mean_sd(.data$Marital_Status == "Separated", denote_sd="paren"),
              "Never Married" = ~ qwraps2::mean_sd(.data$Marital_Status == "Never Married", denote_sd="paren")),
       "Income" = 
         list("Low income" = ~ qwraps2::mean_sd(.data$Income == "Low income", denote_sd="paren"),
              "Middle income" = ~ qwraps2::mean_sd(.data$Income == "Middle income", denote_sd="paren"),
              "High income" = ~ qwraps2::mean_sd(.data$Income == "High income", denote_sd="paren")),
       "Employment Status" = 
         list("Employed for wages" = ~ qwraps2::mean_sd(.data$Employment_Status == "Employed for wages", denote_sd="paren"),
              "Employed for salary" = ~ qwraps2::mean_sd(.data$Employment_Status == "Employed for salary", denote_sd="paren"),
              "Self-employed" = ~ qwraps2::mean_sd(.data$Employment_Status == "Self-employed", denote_sd="paren"),
              "Out of work, looking" = ~ qwraps2::mean_sd(.data$Employment_Status == "Out of work and looking for work", denote_sd="paren"),
              "Out of work, not looking" = ~ qwraps2::mean_sd(.data$Employment_Status == "Out of work but not currently looking for work", denote_sd="paren"),
              "A homemaker" = ~ qwraps2::mean_sd(.data$Employment_Status == "A homemaker", denote_sd="paren"),
              "A student" = ~ qwraps2::mean_sd(.data$Employment_Status == "A student", denote_sd="paren"),
              "Military" = ~ qwraps2::mean_sd(.data$Employment_Status == "Military", denote_sd="paren"),
              "Retired" = ~ qwraps2::mean_sd(.data$Employment_Status == "Retired", denote_sd="paren"),
              "Unable to work" = ~ qwraps2::mean_sd(.data$Employment_Status == "Unable to work", denote_sd="paren"))
  )

bal_tab <- summary_table(dplyr::group_by(ads, Vignette), balance_table)
colnames(bal_tab) <- c("Control", "Bias","Transparency", "Responsiveness",
                 "Control", "Bias","Transparency", "Responsiveness")
covariates <- setNames(c(3, 3, 3, 5, 3, 4, 5, 3, 10), c("Race", "Gender", "Education", "Age", "Pol ID", "AI Knowledge", "Marital Status", "Income", "Employment Status"))
bal_tab <- qable(bal_tab, rgroup=covariates, rnames=rownames(bal_tab), markup="latex", caption="Covariate Balance Table", label="covariate_balance_table") %>% add_header_above(c(" " = 1, "Child Welfare" = 4, "Courts" = 4))
bal_tab #can copy and paste this latex output

save_kable(bal_tab, file="Tables/table_SI1.tex", keep_tex = T)
#only run line above once with final balance table because it throws an irrelevant error that changes the directory and messes up rest of code
```

#####Reproduce Information for Table SI2####
```{r table SI2}
#Create dummy variables for each of the treatments so that we can run regressions
ads <- ads %>% 
          mutate(welfare.control = ifelse(str_detect(Vignette, "Child Welfare Control"), 1, 0),
                 welfare.bias = ifelse(str_detect(Vignette, "Child Welfare Bias"), 1, 0),
                 welfare.transparency = ifelse(str_detect(Vignette, "Child Welfare Transparency"), 1, 0),
                 welfare.responsiveness = ifelse(str_detect(Vignette, "Child Welfare Responsiveness"), 1, 0),
                 courts.control = ifelse(str_detect(Vignette, "Courts Control"), 1, 0),
                 courts.bias = ifelse(str_detect(Vignette, "Courts Bias"), 1, 0),
                 courts.transparency = ifelse(str_detect(Vignette, "Courts Transparency"), 1, 0),
                 courts.responsiveness = ifelse(str_detect(Vignette, "Courts Responsiveness"), 1, 0))
                 
for(level in unique(ads$Vignette)){ #this should do the same as above
  ads[paste(make.names(level))] <- ifelse(ads$Vignette == level, 1, 0)
}  

#Regressions for covariates predicting each of the 8 treatment types
reg.welfare.control <- lm(welfare.control ~ Pol + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data = ads)
reg.welfare.bias <- lm(welfare.bias ~ Pol + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data = ads)
reg.welfare.transparency <- lm(welfare.transparency ~ Pol + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data = ads)
reg.welfare.responsiveness <- lm(welfare.responsiveness ~ Pol + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data = ads)
reg.courts.control <- lm(courts.control ~ Pol + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data = ads)
reg.courts.bias <- lm(courts.bias ~ Pol + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data = ads)
reg.courts.transparency <- lm(courts.transparency ~ Pol + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data = ads)
reg.courts.responsiveness <- lm(courts.responsiveness ~ Pol + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data = ads)

#p-values for global F-tests for each of the 8 models
f.welfare.control <- f_stat_p_value(reg.welfare.control)
f.welfare.bias <- f_stat_p_value(reg.welfare.bias)
f.welfare.transparency <- f_stat_p_value(reg.welfare.transparency)
f.welfare.responsiveness <- f_stat_p_value(reg.welfare.responsiveness)
f.courts.control <- f_stat_p_value(reg.courts.control)
f.courts.bias <- f_stat_p_value(reg.courts.bias)
f.courts.transparency <- f_stat_p_value(reg.courts.transparency)
f.courts.responsiveness <- f_stat_p_value(reg.courts.responsiveness)

f.data.frame <- rbind(f.welfare.control, f.welfare.bias, f.welfare.transparency, f.welfare.responsiveness, f.courts.control, f.courts.bias, f.courts.transparency, f.courts.responsiveness)
row.names(f.data.frame) <- c("Child Welfare Control", "Child Welfare Bias", "Child Welfare Transparency", "Child Welfare Responsiveness", "Courts Control", "Courts Bias", "Courts Transparency", "Courts Responsiveness")
colnames(f.data.frame) <- c("p-value")
f.table <- print(xtable(f.data.frame, type = "latex", digits = 4, caption ="Covariate Balance F-Tests: p-values Across Vignette Treatments", label="covariate_ftest_table"), file="Tables/table_SI2.tex", include.rownames = T, comment=F)
```

#####Reproduce Table SI3####
```{r table SI3}
#Mechanism descriptive outcomes
means <- ads %>% group_by(Mechanism) %>% summarize(Mean_Feeling = mean(Feeling),
                                                         Mean_Trust = mean(Trust),
                                                         Mean_Quality = mean(Quality),
                                                         Mean_Impact = mean(Impact))

sds <- ads %>% group_by(Mechanism) %>% summarize(SD_Feeling = sd(Feeling),
                                                  SD_Trust = sd(Trust),
                                                  SD_Quality = sd(Quality),
                                                  SD_Impact = sd(Impact))

descriptives <- as.data.frame(matrix(rep(0,20), 4, 5))
descriptives[,1] <- means$Mechanism

for(i in 1:4){
  for(j in 2:5){
    descriptives[i,j] <- paste0(round(means[i,j], 2), " (", round(sds[i,j], 2), ")")  
    }
}

colnames(descriptives) <- c("Public Value Failure", "Feeling", "Trust", "Quality", "Impact")

print(descriptives)

print(xtable(descriptives, type="latex", caption="Means and Standard Deviations for Outcomes by Public Value Failure", label="mechanism_descriptives_table"), file="Tables/table_SI3.tex", include.rownames = F, comment=F)
```

#####Reproduce Table SI4####
```{r table SI4}
#Vignette descriptive outcomes
means <- ads %>% group_by(Vignette) %>% summarize(Mean_Feeling = mean(Feeling),
                                                         Mean_Trust = mean(Trust),
                                                         Mean_Quality = mean(Quality),
                                                         Mean_Impact = mean(Impact))

sds <- ads %>% group_by(Vignette) %>% summarize(SD_Feeling = sd(Feeling),
                                                  SD_Trust = sd(Trust),
                                                  SD_Quality = sd(Quality),
                                                  SD_Impact = sd(Impact))

descriptives <- as.data.frame(matrix(rep(0,40), 8, 5))
descriptives[,1] <- c("Child Welfare Control", "Child Welfare Bias", "Child Welfare Transparency", "Child Welfare Responsiveness", "Courts Control", "Courts Bias", "Courts Transparency", "Courts Responsiveness")

for(i in 1:8){
  for(j in 2:5){
    descriptives[i,j] <- paste0(round(means[i,j], 2), " (", round(sds[i,j], 2), ")")  
    }
}

colnames(descriptives) <- c("Vignette", "Feeling", "Trust", "Quality", "Impact")

print(descriptives)

print(xtable(descriptives, type="latex", caption="Means and Standard Deviations for Outcomes by Vignette", label="vignette_descriptives_table"), file="Tables/table_SI4.tex", include.rownames = F, comment=F)
```

####Reproduce Figure SI3####
```{r figure SI3}
#Create figure
vignette_outcomes <- ads %>% group_by(Vignette) %>% summarize(
  mean_feeling = mean(Feeling), sd_feeling = sd(Feeling),
  mean_trust = mean(Trust), sd_trust = sd(Trust),
  mean_quality = mean(Quality), sd_quality = sd(Quality),
  mean_impact = mean(Impact), sd_impact = sd(Impact))

vignette_outcomes <- vignette_outcomes %>% pivot_longer(cols = 2:9) %>% 
  mutate(outcome = rep(rep(c("Feeling", "Trust", "Quality", "Impact"), each = 2),8)) %>% 
  ungroup %>% 
  mutate(name = rep(c("mean","sd"),32)) %>% 
  pivot_wider(id_cols = c(1,4), names_from = name, values_from = value)
vignette_outcomes$outcome <- factor(vignette_outcomes$outcome, levels = c("Feeling", "Trust", "Quality", "Impact"))


feeling_vignette <- ggplot(vignette_outcomes %>% filter(outcome == "Feeling"), 
                      aes(x=fct_rev(Vignette), y=mean)) + 
  geom_point(size = 1.5) +
  geom_hline(yintercept=0, colour = "grey60", lty=2) + 
  coord_flip() +
  xlab("") + 
  ylab("\n") +
  labs(title="Feeling") +
  ylim(-5,5) + 
  theme_bw() + 
  theme(plot.title = element_text(face="bold", hjust = 0.5),
        #legend.position = "none",
        legend.justification=c(.6, 0),
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank()
  )

trust_vignette <- ggplot(vignette_outcomes %>% filter(outcome == "Trust"), 
                    aes(x=fct_rev(Vignette), y=mean)) + 
  geom_point(size = 1.5) +
  geom_hline(yintercept=0, colour = "grey60", lty=2) + 
  coord_flip() +
  xlab("") + 
  ylab("\n") +
  labs(title="Trust") +
  ylim(-5,5) + 
  theme_bw() + 
  theme(plot.title = element_text(face="bold", hjust = 0.5),
        legend.position = "none",
        #legend.justification=c(.6, 0),
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank()
  )

quality_vignette <- ggplot(vignette_outcomes %>% filter(outcome == "Quality"), 
                      aes(x=fct_rev(Vignette), y=mean)) + 
  geom_point(size = 1.5) +
  geom_hline(yintercept=0, colour = "grey60", lty=2) + 
  coord_flip() +
  xlab("") + 
  ylab("\nMean Outcome") +
  labs(title="Quality") +
  ylim(-5,5) + 
  theme_bw() + 
  theme(plot.title = element_text(face="bold", hjust = 0.5),
        legend.position = "none",
        #legend.justification=c(.6, 0),
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank()
  )

impact_vignette <- ggplot(vignette_outcomes %>% filter(outcome == "Impact"), 
                     aes(x=fct_rev(Vignette), y=mean,)) + 
  geom_point(size = 1.5) +
  geom_hline(yintercept=4, colour = "grey60", lty=2) + 
  coord_flip() +
  xlab("") + 
  ylab("\nMean Outcome") +
  labs(title="Impact") +
  ylim(2.5,5.5) + 
  theme_bw() + 
  theme(plot.title = element_text(face="bold", hjust = 0.5),
        legend.position = "none",
        #legend.justification=c(.6, 0),
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank()
  )

jpeg("Figures/fig_SI3.jpg", units="in", width=9, height=7, res=300)
ggarrange(feeling_vignette, trust_vignette, quality_vignette, impact_vignette, nrow=2, ncol = 2, common.legend = T, legend = "bottom")
dev.off()
```

####Reproduce Table SI5####
```{r table SI5}
out1_unadjust <- lm(Feeling ~ Mechanism, data=ads)
out2_unadjust <- lm(Trust ~ Mechanism, data=ads)
out3_unadjust <- lm(Quality ~ Mechanism, data=ads)
out4_unadjust <- lm(Impact ~ Mechanism, data=ads)

stargazer(out1_unadjust, out2_unadjust, out3_unadjust, out4_unadjust,
          se=list(sqrt(diag(vcovHC(out1_unadjust,type="HC1"))), 
                  sqrt(diag(vcovHC(out2_unadjust,type="HC1"))),
                  sqrt(diag(vcovHC(out3_unadjust,type="HC1"))),
                  sqrt(diag(vcovHC(out4_unadjust,type="HC1")))),
          covariate.labels=c("Bias", "Lack of Transparency", "Lack of Responsiveness"),
          title="Public Value Failure Regression Results (Unadjusted)",
          label="mechanism_table_unadjusted",
          notes="With robust SEs, excluding covariates, and without block coefficients printed",
          style="APSR",
          header=F,
          type="latex",
          out="Tables/table_SI5.tex")
```

####Reproduce Information for Table SI6####
```{r table SI6}
#First, does the probability of treatment vary by block?  This is helpful to assess the potential for bias.
block_summary <- ads %>% group_by(Block_ID) %>%
  summarize(n = length(Block_ID),
            prop = length(Block_ID)/nrow(ads),
            control_prop = length(which(Mechanism=="Control"))/length(Block_ID),
            bias_prop = length(which(Mechanism=="Bias"))/length(Block_ID),
            trans_prop = length(which(Mechanism=="Transparency"))/length(Block_ID),
            resp_prop = length(which(Mechanism=="Responsiveness"))/length(Block_ID)
  )
#Probability of assignment to each of the main treatment groups is pretty much the same across groups,
#except for groups that are very small (note: this was of course by design!)

#Now calculate block-weighted ATEs
block_ATEs <- ads %>% group_by(Block_ID, Mechanism) %>%
  summarize(feelings_avg = mean(Feeling, na.rm=T),
            trust_avg = mean(Trust, na.rm=T),
            quality_avg = mean(Quality, na.rm=T),
            impact_avg = mean(Impact, na.rm=T)
            )
block_controls <- block_ATEs %>% filter(Mechanism=="Control") %>% pivot_longer(3:6)
block_treatments <- block_ATEs %>% filter(Mechanism!="Control") %>% pivot_longer(3:6)
block_ATEs <- left_join(block_treatments, block_controls, by = c("Block_ID","name")) %>%
  mutate(block_ate = value.x - value.y) %>%
  left_join(block_summary[,c(1,3)], by="Block_ID")
block_ATEs <- block_ATEs %>% ungroup() %>% 
  mutate(Outcome = rep(c("Feelings", "Trust", "Quality", "Impact"),56)) %>%
  select(Block_ID, prop, Outcome, Mechanism = Mechanism.x, block_ate)

block_weighted_ates <- block_ATEs %>% group_by(Mechanism, Outcome) %>%
  summarize(block_weighted_ate = weighted.mean(block_ate, prop, na.rm=T))

#Now calculate SEs
block_vars <- ads %>% group_by(Block_ID, Mechanism) %>%
  summarize(feelings_var = var(Feeling, na.rm=T),
            trust_var = var(Trust, na.rm=T),
            quality_var = var(Quality, na.rm=T),
            impact_var = var(Impact, na.rm=T),
            n = n()
            )
block_controls <- block_vars %>% filter(Mechanism=="Control") %>% pivot_longer(3:6)
block_treatments <- block_vars %>% filter(Mechanism!="Control") %>% pivot_longer(3:6)
block_vars <- left_join(block_treatments, block_controls, by = c("Block_ID","name")) %>%
  mutate(block_se =  sqrt((value.x/n.x) + (value.y/n.y))) %>%
  left_join(block_summary[,c(1,3)], by="Block_ID")
block_ses <- block_vars %>% ungroup() %>% 
  mutate(Outcome = rep(c("Feelings", "Trust", "Quality", "Impact"),56)) %>%
  select(Block_ID, prop, Outcome, Mechanism = Mechanism.x, block_se)

block_weighted_ses <- block_ses %>% mutate(part = (block_se^2)*(prop^2)) %>%
  group_by(Mechanism, Outcome) %>%
  summarize(block_weighted_se = sqrt(sum(part, na.rm=T)))

block_values <- block_weighted_ates %>% ungroup %>% mutate(block_weighted_se = block_weighted_ses[,3] %>% pull)

#Information for Table SI6
block_values %>% mutate(zscore = block_weighted_ate/block_weighted_se, pval = 2*(1-pnorm(abs(zscore), lower.tail = T)))
```

####Reproduce Information for Table SI7####
```{r table SI7}
#Perform regressions
out1_adjust <- lm(Feeling ~ Vignette + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)
out2_adjust <- lm(Trust ~ Vignette + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)
out3_adjust <- lm(Quality ~ Vignette + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)
out4_adjust <- lm(Impact ~ Vignette + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)

#Create table
stargazer(out1_adjust, out2_adjust, out3_adjust, out4_adjust,
          se=list(sqrt(diag(vcovHC(out1_adjust,type="HC1"))), 
                  sqrt(diag(vcovHC(out2_adjust,type="HC1"))),
                  sqrt(diag(vcovHC(out3_adjust,type="HC1"))),
                  sqrt(diag(vcovHC(out4_adjust,type="HC1")))),
          omit=c("GenderOther", "Marital", "Employment", "Industry"),
          covariate.labels=c("Child Welfare Bias", "Child Welfare Lack of Transparency", "Child Welfare Lack of Responsiveness",
                             "Courts Control", "Courts Bias", "Courts Lack of Transparency", 
                             "Courts Lack of Responsiveness",
                             "Liberal", "Conservative", "Low AI Knowledge", "Some AI Knowledge", "High AI Knowledge",
                             "Female", "Millenial", "Gen X", "Boomer", "Silent Gen", "Black", "Other Race", "Some College",
                             "Bachelor's or Higher", "Middle Income", "High Income"),
          title="Vignette Regression Results (Covariate-Adjusted)",
          label="vignette_table_adjusted",
          notes="With robust SEs, including covariates, and without block coefficients and some covariates printed",
          style="APSR",
          header=F,
          type="latex",
          out="Tables/table_SI7.tex")
```

####Reproduce Table SI8####
```{r table SI8}
#Perform regressions
out1_unadjust <- lm(Feeling ~ Vignette, data=ads)
out2_unadjust <- lm(Trust ~ Vignette, data=ads)
out3_unadjust <- lm(Quality ~ Vignette, data=ads)
out4_unadjust <- lm(Impact ~ Vignette, data=ads)

#Create table
stargazer(out1_unadjust, out2_unadjust, out3_unadjust, out4_unadjust,
          se=list(sqrt(diag(vcovHC(out1_unadjust,type="HC1"))), 
                  sqrt(diag(vcovHC(out2_unadjust,type="HC1"))),
                  sqrt(diag(vcovHC(out3_unadjust,type="HC1"))),
                  sqrt(diag(vcovHC(out4_unadjust,type="HC1")))),
          covariate.labels=c("Child Welfare Bias", "Child Welfare Lack of Transparency", "Child Welfare Lack of Responsiveness",
                             "Courts Control", "Courts Bias", "Courts Lack of Transparency", 
                             "Courts Lack of Responsiveness"),
          title="Vignette Regression Results (Unadjusted)",
          label="vignette_table_unadjusted",
          notes="With robust SEs, excluding covariates, and without block coefficients printed",
          style="APSR",
          header=F,
          type="latex",
          out="Tables/table_SI8.tex")
```

####Reproduce Table SI9####
```{r table SI9}
#Perform regressions
out1_unadjust_conservative <- lm(Feeling ~ Vignette, data=ads %>% filter(Pol=="Conservative"))
out2_unadjust_conservative <- lm(Trust ~ Vignette, data=ads %>% filter(Pol=="Conservative"))
out3_unadjust_conservative <- lm(Quality ~ Vignette, data=ads %>% filter(Pol=="Conservative"))
out4_unadjust_conservative <- lm(Impact ~ Vignette, data=ads %>% filter(Pol=="Conservative"))

#Create table
stargazer(out1_unadjust_conservative, out2_unadjust_conservative, out3_unadjust_conservative, out4_unadjust_conservative,
           se=list(sqrt(diag(vcovHC(out1_unadjust_conservative,type="HC1"))), 
                  sqrt(diag(vcovHC(out2_unadjust_conservative,type="HC1"))),
                  sqrt(diag(vcovHC(out3_unadjust_conservative,type="HC1"))),
                  sqrt(diag(vcovHC(out4_unadjust_conservative,type="HC1")))),
          covariate.labels=c("Child Welfare Bias", "Child Welfare Lack of Transparency", "Child Welfare Lack of Responsiveness",
                             "Courts Control", "Courts Bias", "Courts Lack of Transparency", 
                             "Courts Lack of Responsiveness"),
          title="Conservative Subgroup Regression Results (Unadjusted)",
          label="conservative_table_unadjusted",
          notes="With robust SEs, excluding covariates, and without block coefficients printed",
          style="APSR",
          header=F,
          type="latex",
          out="Tables/table_SI9.tex")
```

####Reproduce Table SI10####
```{r table SI10}
#Perform regressions
out1_unadjust_liberal <- lm(Feeling ~ Vignette, data=ads %>% filter(Pol=="Liberal"))
out2_unadjust_liberal <- lm(Trust ~ Vignette, data=ads %>% filter(Pol=="Liberal"))
out3_unadjust_liberal <- lm(Quality ~ Vignette, data=ads %>% filter(Pol=="Liberal"))
out4_unadjust_liberal <- lm(Impact ~ Vignette, data=ads %>% filter(Pol=="Liberal"))

#Create table
stargazer(out1_unadjust_liberal, out2_unadjust_liberal, out3_unadjust_liberal, out4_unadjust_liberal,
           se=list(sqrt(diag(vcovHC(out1_unadjust_liberal,type="HC1"))), 
                  sqrt(diag(vcovHC(out2_unadjust_liberal,type="HC1"))),
                  sqrt(diag(vcovHC(out3_unadjust_liberal,type="HC1"))),
                  sqrt(diag(vcovHC(out4_unadjust_liberal,type="HC1")))),
          covariate.labels=c("Child Welfare Bias", "Child Welfare Lack of Transparency", "Child Welfare Lack of Responsiveness",
                             "Courts Control", "Courts Bias", "Courts Lack of Transparency", 
                             "Courts Lack of Responsiveness"),
          title="Liberal Subgroup Regression Results (Unadjusted)",
          label="liberal_table_unadjusted",
          notes="With robust SEs, excluding covariates, and without block coefficients printed",
          style="APSR",
          header=F,
          type="latex",
          out="Tables/table_SI10.tex")
```

#Whether to Standardize Outcome Variables
```{r standardize variables}
standardize <- TRUE

#because SDs and sample sizes are similar across treatment groups, we use Cohen's d (scale outcomes by sd of whole sample)
if(standardize == TRUE) {
  ads$Feeling <- scale(ads$Feeling)
  ads$Trust <- scale(ads$Trust)
  ads$Quality <- scale(ads$Quality)
  ads$Impact <- scale(ads$Impact)
  ads$Petition <- scale(ads$Petition)
  ads$Meeting <- scale(ads$Meeting)
}
```

####Reproduce Figure SI4####
```{r figure SI4}
#Perform regressions for conservatives
out1_adjust_conservative <- lm(Feeling ~ Vignette + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Conservative"))
out2_adjust_conservative <- lm(Trust ~ Vignette + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Conservative")) 
out3_adjust_conservative <- lm(Quality ~ Vignette + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Conservative"))
out4_adjust_conservative <- lm(Impact ~ Vignette + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Conservative"))

out1_adjust_conservative_tidy <- tidy(out1_adjust_conservative)[2:4,] %>% mutate(model = "Child Welfare")
out2_adjust_conservative_tidy <- tidy(out2_adjust_conservative)[2:4,] %>% mutate(model = "Child Welfare")
out3_adjust_conservative_tidy <- tidy(out3_adjust_conservative)[2:4,] %>% mutate(model = "Child Welfare")
out4_adjust_conservative_tidy <- tidy(out4_adjust_conservative)[2:4,] %>% mutate(model = "Child Welfare")

out1_adjust_conservative_opp <- lm(Feeling ~ relevel(Vignette, ref="Courts Control") + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Conservative"))
out1_adjust_liberal_opp <- lm(Feeling ~ relevel(Vignette, ref="Courts Control") + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Liberal"))
out2_adjust_conservative_opp <- lm(Trust ~ relevel(Vignette, ref="Courts Control") + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Conservative")) 
out2_adjust_liberal_opp <- lm(Trust ~ relevel(Vignette, ref="Courts Control") + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Liberal"))
out3_adjust_conservative_opp <- lm(Quality ~ relevel(Vignette, ref="Courts Control") + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Conservative"))
out3_adjust_liberal_opp <- lm(Quality ~ relevel(Vignette, ref="Courts Control") + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Liberal"))
out4_adjust_conservative_opp <- lm(Impact ~ relevel(Vignette, ref="Courts Control") + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Conservative"))
out4_adjust_liberal_opp <- lm(Impact ~ relevel(Vignette, ref="Courts Control") + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Liberal"))

names(out1_adjust_conservative_opp$coefficients)[6:8] <- c("VignetteCourts Bias", "VignetteCourts Transparency", "VignetteCourts Responsiveness")
names(out2_adjust_conservative_opp$coefficients)[6:8] <- c("VignetteCourts Bias", "VignetteCourts Transparency", "VignetteCourts Responsiveness")
names(out3_adjust_conservative_opp$coefficients)[6:8] <- c("VignetteCourts Bias", "VignetteCourts Transparency", "VignetteCourts Responsiveness")
names(out4_adjust_conservative_opp$coefficients)[6:8] <- c("VignetteCourts Bias", "VignetteCourts Transparency", "VignetteCourts Responsiveness")

out1_adjust_conservative_opp_tidy <- tidy(out1_adjust_conservative_opp)[6:8,] %>% mutate(model = "Courts")
out2_adjust_conservative_opp_tidy <- tidy(out2_adjust_conservative_opp)[6:8,] %>% mutate(model = "Courts")
out3_adjust_conservative_opp_tidy <- tidy(out3_adjust_conservative_opp)[6:8,] %>% mutate(model = "Courts")
out4_adjust_conservative_opp_tidy <- tidy(out4_adjust_conservative_opp)[6:8,] %>% mutate(model = "Courts")

feeling_con <- rbind(out1_adjust_conservative_tidy, out1_adjust_conservative_opp_tidy) %>% arrange((model))
feeling_con$term <- rep(c("Bias","Lack of Transparency","Lack of Responsiveness"),2)

trust_con <- rbind(out2_adjust_conservative_tidy, out2_adjust_conservative_opp_tidy)
trust_con$term <- rep(c("Bias","Lack of Transparency","Lack of Responsiveness"),2)

quality_con <- rbind(out3_adjust_conservative_tidy, out3_adjust_conservative_opp_tidy)
quality_con$term <- rep(c("Bias","Lack of Transparency","Lack of Responsiveness"),2)

impact_con <- rbind(out4_adjust_conservative_tidy, out4_adjust_conservative_opp_tidy)
impact_con$term <- rep(c("Bias","Lack of Transparency","Lack of Responsiveness"),2)

#Create figure
jpeg("Figures/fig_SI4.jpg", units="in", width=9, height=7, res=300)

con_feeling <- dwplot(arrange(feeling_con, desc(model)),
                      vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),
                      dot_args = list(aes(shape = model)), size=1.5) +
  scale_color_manual(name = "model", breaks=c("Child Welfare","Courts"), values = c("#C77CFF","#7CAE00")) +
  scale_shape_discrete(name = "model", breaks=c("Child Welfare","Courts")) +
  xlab("\n ") + ylab("") +
  xlim(-1.0, .5) +
  ggtitle("Feeling") +
  theme_bw() + 
  theme(plot.title = element_text(face="bold", hjust = 0.5),
        legend.justification=c(.6, 0),
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank()
  )

con_trust <- dwplot(arrange(trust_con, desc(model)),
                    vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),
                    dot_args = list(aes(shape = model)), size=1.5) +
  scale_color_manual(name = "model", breaks=c("Child Welfare","Courts"), values = c("#C77CFF","#7CAE00")) +
  scale_shape_discrete(name = "model", breaks=c("Child Welfare","Courts")) +
  xlab("\n ") + ylab("") +
  xlim(-1.0, .5) +
  ggtitle("Trust") +
  theme_bw() + 
  theme(plot.title = element_text(face="bold", hjust = 0.5),
        legend.position = "none",
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank())

con_quality <- dwplot(arrange(quality_con, desc(model)),
                      vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),
                      dot_args = list(aes(shape = model)), size=1.5) +
  scale_color_manual(name = "model", breaks=c("Child Welfare","Courts"), values = c("#C77CFF","#7CAE00")) +
  scale_shape_discrete(name = "model", breaks=c("Child Welfare","Courts")) +
  xlab("\nStandardized Treatment Effect") + ylab("") +
  xlim(-1.0, .5) +
  ggtitle("Quality") +
  theme_bw() + 
  theme(plot.title = element_text(face="bold", hjust = 0.5),
        legend.position = "none",
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank())

con_impact <- dwplot(arrange(impact_con, desc(model)),
                     vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),
                     dot_args = list(aes(shape = model)), size=1.5) +
  scale_color_manual(name = "model", breaks=c("Child Welfare","Courts"), values = c("#C77CFF","#7CAE00")) +
  scale_shape_discrete(name = "model", breaks=c("Child Welfare","Courts")) +
  xlab("\nStandardized Treatment Effect") + ylab("") +
  xlim(-1.0, .5) +
  ggtitle("Impact") +
  theme_bw() + 
  theme(plot.title = element_text(face="bold", hjust = 0.5),
        legend.position = "none",
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank())

ggarrange(con_feeling, con_trust, con_quality, con_impact, nrow=2, ncol = 2, common.legend = T, legend = "bottom")

dev.off()
```

####Reproduce Figure SI5####
```{r figure SI5}
#Perform regressions for liberals
out1_adjust_liberal <- lm(Feeling ~ Vignette + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Liberal"))
out2_adjust_liberal <- lm(Trust ~ Vignette + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Liberal"))
out3_adjust_liberal <- lm(Quality ~ Vignette + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Liberal"))
out4_adjust_liberal <- lm(Impact ~ Vignette + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Liberal"))

out1_adjust_liberal_tidy <- tidy(out1_adjust_liberal)[2:4,] %>% mutate(model = "Child Welfare")
out2_adjust_liberal_tidy <- tidy(out2_adjust_liberal)[2:4,] %>% mutate(model = "Child Welfare")
out3_adjust_liberal_tidy <- tidy(out3_adjust_liberal)[2:4,] %>% mutate(model = "Child Welfare")
out4_adjust_liberal_tidy <- tidy(out4_adjust_liberal)[2:4,] %>% mutate(model = "Child Welfare")

out1_adjust_liberal_opp <- lm(Feeling ~ relevel(Vignette, ref="Courts Control") + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Liberal"))
out2_adjust_liberal_opp <- lm(Trust ~ relevel(Vignette, ref="Courts Control") + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Liberal"))
out3_adjust_liberal_opp <- lm(Quality ~ relevel(Vignette, ref="Courts Control") + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Liberal"))
out4_adjust_liberal_opp <- lm(Impact ~ relevel(Vignette, ref="Courts Control") + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Liberal"))

out1_adjust_liberal_opp_tidy <- tidy(out1_adjust_liberal_opp)[6:8,] %>% mutate(model = "Courts")
out2_adjust_liberal_opp_tidy <- tidy(out2_adjust_liberal_opp)[6:8,] %>% mutate(model = "Courts")
out3_adjust_liberal_opp_tidy <- tidy(out3_adjust_liberal_opp)[6:8,] %>% mutate(model = "Courts")
out4_adjust_liberal_opp_tidy <- tidy(out4_adjust_liberal_opp)[6:8,] %>% mutate(model = "Courts")

feeling_lib <- rbind(out1_adjust_liberal_tidy, out1_adjust_liberal_opp_tidy) %>% arrange((model))
feeling_lib$term <- rep(c("Bias","Lack of Transparency","Lack of Responsiveness"),2)

trust_lib <- rbind(out2_adjust_liberal_tidy, out2_adjust_liberal_opp_tidy)
trust_lib$term <- rep(c("Bias","Lack of Transparency","Lack of Responsiveness"),2)

quality_lib <- rbind(out3_adjust_liberal_tidy, out3_adjust_liberal_opp_tidy)
quality_lib$term <- rep(c("Bias","Lack of Transparency","Lack of Responsiveness"),2)

impact_lib <- rbind(out4_adjust_liberal_tidy, out4_adjust_liberal_opp_tidy)
impact_lib$term <- rep(c("Bias","Lack of Transparency","Lack of Responsiveness"),2)

names(out1_adjust_liberal_opp$coefficients)[6:8] <- c("VignetteCourts Bias", "VignetteCourts Transparency", "VignetteCourts Responsiveness")
names(out2_adjust_liberal_opp$coefficients)[6:8] <- c("VignetteCourts Bias", "VignetteCourts Transparency", "VignetteCourts Responsiveness")
names(out3_adjust_liberal_opp$coefficients)[6:8] <- c("VignetteCourts Bias", "VignetteCourts Transparency", "VignetteCourts Responsiveness")
names(out4_adjust_liberal_opp$coefficients)[6:8] <- c("VignetteCourts Bias", "VignetteCourts Transparency", "VignetteCourts Responsiveness")

#Create figure
jpeg("Figures/fig_SI5.jpg", units="in", width=9, height=7, res=300)

lib_feeling <- dwplot(arrange(feeling_lib, desc(model)),
                      vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),
                      dot_args = list(aes(shape = model)), size=1.5) +
  scale_color_manual(name = "model", breaks=c("Child Welfare","Courts"), values = c("#C77CFF","#7CAE00")) +
  scale_shape_discrete(name = "model", breaks=c("Child Welfare","Courts")) +
  xlab("\n ") + ylab("") +
  xlim(-1.0, .5) +
  ggtitle("Feeling") +
  theme_bw() + 
  theme(plot.title = element_text(face="bold", hjust = 0.5),
        legend.justification=c(.6, 0),
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank()
  )

lib_trust <- dwplot(arrange(trust_lib, desc(model)),
                    vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),
                    dot_args = list(aes(shape = model)), size=1.5) +
  scale_color_manual(name = "model", breaks=c("Child Welfare","Courts"), values = c("#C77CFF","#7CAE00")) +
  scale_shape_discrete(name = "model", breaks=c("Child Welfare","Courts")) +
  xlab("\n ") + ylab("") +
  xlim(-1.0, .5) +
  ggtitle("Trust") +
  theme_bw() + 
  theme(plot.title = element_text(face="bold", hjust = 0.5),
        legend.position = "none",
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank())

lib_quality <- dwplot(arrange(quality_lib, desc(model)),
                      vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),
                      dot_args = list(aes(shape = model)), size=1.5) +
  scale_color_manual(name = "model", breaks=c("Child Welfare","Courts"), values = c("#C77CFF","#7CAE00")) +
  scale_shape_discrete(name = "model", breaks=c("Child Welfare","Courts")) +
  xlab("\nStandardized Treatment Effect") + ylab("") +
  xlim(-1.0, .5) +
  ggtitle("Quality") +
  theme_bw() + 
  theme(plot.title = element_text(face="bold", hjust = 0.5),
        legend.position = "none",
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank())

lib_impact <- dwplot(arrange(impact_lib, desc(model)),
                     vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),
                     dot_args = list(aes(shape = model)), size=1.5) +
  scale_color_manual(name = "model", breaks=c("Child Welfare","Courts"), values = c("#C77CFF","#7CAE00")) +
  scale_shape_discrete(name = "model", breaks=c("Child Welfare","Courts")) +
  xlab("\nStandardized Treatment Effect") + ylab("") +
  xlim(-1.0, .5) +
  ggtitle("Impact") +
  theme_bw() + 
  theme(plot.title = element_text(face="bold", hjust = 0.5),
        legend.position = "none",
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank())

ggarrange(lib_feeling, lib_trust, lib_quality, lib_impact, nrow=2, ncol = 2, common.legend = T, legend = "bottom")

dev.off()
```

####Reproduce Figure SI6####
```{r figure SI6}
#Perform regressions
out1_adjust <- lm(Petition ~ Mechanism + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)
out2_adjust <- lm(Meeting ~ Mechanism + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)

out1_adjust_tidy <- tidy(out1_adjust)[2:4,] %>% mutate(model = "Petition")
out2_adjust_tidy <- tidy(out2_adjust)[2:4,] %>% mutate(model = "Meeting")

main <- rbind(out1_adjust_tidy, out2_adjust_tidy) %>% arrange((model))
main$term <- rep(c("Bias","Lack of Transparency","Lack of Responsiveness"),2)
main$model <- factor(main$model, levels = c("Petition","Meeting"))
main <- main %>% arrange(model)

#Create figure
jpeg("Figures/fig_SI6.jpg", units="in", width=9, height=7, res=300)

exploratory_main <- dwplot(main %>% arrange(desc(model)),
vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),
dot_args = list(aes(shape = model), size = 1.5)) +
scale_color_discrete(name = "Outcome", breaks=c("Petition","Meeting")) +
  scale_shape_discrete(name = "Outcome", breaks=c("Petition","Meeting")) +
    xlab("\nStandardized Treatment Effect") + ylab("") +
    xlim(-.75, .75) +
    theme_bw() + 
    theme(plot.title = element_text(face="bold", hjust = 0.5),
          legend.position = "right",
          legend.background = element_rect(colour="grey80")
          )

exploratory_main
dev.off()
```

####Reproduce Figure SI7####
```{r figure SI7}
#Perform regressions
out1_adjust <- lm(Petition ~ Vignette + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)
out2_adjust <- lm(Meeting ~ Vignette + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)

out1_adjust_tidy <- tidy(out1_adjust)[2:4,] %>% mutate(model = "Child Welfare")
out2_adjust_tidy <- tidy(out2_adjust)[2:4,] %>% mutate(model = "Child Welfare")

out1_adjust_opp <- lm(Petition ~ relevel(Vignette, ref="Courts Control") + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)
out2_adjust_opp <- lm(Meeting ~ relevel(Vignette, ref="Courts Control") + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)
names(out1_adjust_opp$coefficients)[6:8] <- c("VignetteCourts Bias", "VignetteCourts Transparency", "VignetteCourts Responsiveness")
names(out2_adjust_opp$coefficients)[6:8] <- c("VignetteCourts Bias", "VignetteCourts Transparency", "VignetteCourts Responsiveness")

out1_adjust_opp_tidy <- tidy(out1_adjust_opp)[6:8,] %>% mutate(model = "Courts")
out2_adjust_opp_tidy <- tidy(out2_adjust_opp)[6:8,] %>% mutate(model = "Courts")

petition <- rbind(out1_adjust_tidy, out1_adjust_opp_tidy) %>% arrange((model))
petition$term <- rep(c("Bias","Lack of Transparency","Lack of Responsiveness"),2)

meeting <- rbind(out2_adjust_tidy, out2_adjust_opp_tidy)
meeting$term <- rep(c("Bias","Lack of Transparency","Lack of Responsiveness"),2)

#Create figure
jpeg("Figures/fig_SI7.jpg", units="in", width=9, height=3.5, res=300)

dw_petition <- dwplot(arrange(petition, desc(model)),
                     vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),
                     dot_args = list(aes(shape = model)), size=1.5) +
  scale_color_manual(name = "model", breaks=c("Child Welfare","Courts"), values = c("#C77CFF","#7CAE00")) +
  scale_shape_discrete(name = "model", breaks=c("Child Welfare","Courts")) +
  xlab("\nStandardized Treatment Effect") + ylab("") +
  xlim(-.75, .75) +
  ggtitle("Petition") +
  theme_bw() + 
  theme(plot.title = element_text(face="bold", hjust = 0.5),
        legend.justification=c(.6, 0),
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank()
  )

dw_meeting <- dwplot(arrange(meeting, desc(model)),
vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),
                     dot_args = list(aes(shape = model)), size=1.5) +
  scale_color_manual(name = "model", breaks=c("Child Welfare","Courts"), values = c("#C77CFF","#7CAE00")) +
  scale_shape_discrete(name = "model", breaks=c("Child Welfare","Courts")) +
    xlab("\nStandardized Treatment Effect") + ylab("") +
    xlim(-.75, .75) +
    ggtitle("Meeting") +
    theme_bw() + 
    theme(plot.title = element_text(face="bold", hjust = 0.5),
          legend.position = "none",
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank())

ggarrange(dw_petition, dw_meeting, nrow=1, ncol = 2, common.legend = T, legend = "bottom")

dev.off()
```

#####Reproduce Information for Table SI11####
```{r table SI11}
load("Data/ads_clean_attentive.RData")

#Perform regressions
out1_adjust <- lm(Feeling ~ Mechanism + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)
out2_adjust <- lm(Trust ~ Mechanism + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)
out3_adjust <- lm(Quality ~ Mechanism + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)
out4_adjust <- lm(Impact ~ Mechanism + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)

#Create table
stargazer(out1_adjust, out2_adjust, out3_adjust, out4_adjust,
          se=list(sqrt(diag(vcovHC(out1_adjust,type="HC1"))), 
                  sqrt(diag(vcovHC(out2_adjust,type="HC1"))),
                  sqrt(diag(vcovHC(out3_adjust,type="HC1"))),
                  sqrt(diag(vcovHC(out4_adjust,type="HC1")))),
          omit=c("GenderOther", "Marital", "Employment", "Industry"),
          covariate.labels=c("Bias", "Lack of Transparency", "Lack of Responsiveness",
                             "Liberal", "Conservative", "Low AI Knowledge", "Some AI Knowledge", "High AI Knowledge",
                             "Female", "Millenial", "Gen X", "Boomer", "Silent Gen", "Black", "Other Race", "Some College",
                             "Bachelor's or Higher", "Middle Income", "High Income"),
          title="Public Value Failure Regression Results (Covariate-Adjusted)",
          label="mechanism_table_adjusted",
          notes="With robust SEs, including covariates, and without block coefficients and some covariates printed",
          style="APSR",
          header=F,
          type="latex",
          out="Tables/table_SI11.tex")
```

